//#include <bits/stl_alloc.h>
#include <stdexcept>
#include <string>
#include <algorithm>
#include <iostream>
#include <utility>
#include <vector>
#include <sstream>
#include <fstream>
#include <deque>
#include <boost/lexical_cast.hpp>
#include <boost/tuple/tuple.hpp>

# define DEBUG_LEVEL 1
#include <debug_func.h>
#include <perlutil.h>
#include <ultoa.h>
#include <bl2seq_mine.h>
#include <portable_timer.h>
#include <progress_indicator.h>
#include <bcpp_parsefasta.h>
#include <bcpp_read_seq.h>
#include <stlhelper.h>
#include <seq_pairs.h>
#include <string_hash.h>
#include <print_error.h>
#include <bcpp_read_seq.h>

#include "get_arg.h"

using BioCPP::t_collect_sequences;
using BioCPP::t_collect_sequences_report_errors_to_strm;
using BioCPP::eSEQ_NODUPLICATES;
using BioCPP::eSEQ_DUPLICATES;
using std::cerr;
using std::cout;
using std::string;
using std::vector;
using std::deque;
using std::istream;
using std::ifstream;
using std::sort;
using std::make_pair;
using std::pair;
using std::ostringstream;
using boost::lexical_cast;
using boost::tie;
using bcpp_align::t_align_results;
//________________________________________________________________________________________

//		data structures, hashes etc.

//________________________________________________________________________________________

struct t_align_error
{
	string name1;
	string name2;
	string error_str;
	t_align_error(const string& n1, const string& n2, const string& err)
			:name1(n1), name2(n2), error_str(err){}
};
inline std::ostream& operator << (std::ostream& ostrm, const t_align_error& errors)
{
	ostrm
		<< errors.name1                         << "\t"
		<< errors.name2                         << "\t"
		<< errors.error_str						<< "\n";
	return ostrm;
}

struct t_align_results_pair
{
	string name1;
	string name2;
	t_align_results blast_result;
	t_align_results_pair(const string& n1, const string& n2, const t_align_results& blast_result_)
			:name1(n1), name2(n2), blast_result(blast_result_){}
};

inline std::ostream& operator << (std::ostream& ostrm, const t_align_results_pair& results)
{
	ostrm
		<< results.name1                            << "\t"
		<< results.name2                            << "\t"
		<< results.blast_result						<< "\n";
	return ostrm;
}


//________________________________________________________________________________________

//		align_sequences

//________________________________________________________________________________________
bool align_sequences(	t_progress_indicator&	dots,
						t_blast_options& 		blast_options,
						unsigned				lencutoff,
						t_collect_sequences&	seq_db,
						unsigned				seq1,
						unsigned				seq2,
						t_program_args&			args)
{

	++dots;

	debug_func_cerr( 8 );

	t_align_results curr_results;
	string name1 = seq_db.name(seq1);
	string name2 = seq_db.name(seq2);


	try
	{
		if (do_bl2seq(	blast_options,
						seq_db.sequence(seq1),
						seq_db.sequence(seq2),
						curr_results))
		{
			if (curr_results.align_len < lencutoff)
			{
				args.BLAST_LOG()	<< name1 << "\t"
									<< name2 << "\t"
									<< curr_results.error_message
									<< "Alignment too short. (sequence lengths = "
									<< curr_results.len1 << " and "
									<< curr_results.len2 << ". Aligned positions = "
									<< curr_results.align_len << ", residues = "
									<< curr_results.aligned_residues << ")\n";
				return false;
			}
			else
			{
				cout << t_align_results_pair(
									name1,
									name2,
									curr_results);
				return true;
			}
		}
	}
	catch (std::exception& e)
	{
		curr_results.error_message += e.what();
	};
	
	args.BLAST_LOG() << (t_align_error(
						name1,
						name2,
						"No similarity. " +
						curr_results.error_message));
	return false;
}

















//________________________________________________________________________________________

//		read_seqid_identifiers

//			for each identifier, construct all vs this identifier

//________________________________________________________________________________________
bool read_seqid_identifiers(istream& 							f_seq_names,
							const t_collect_sequences&			seq_db,
							vector<pair<unsigned, unsigned> >& 	sequence_pairs,
							ostream&							CERR,
							bool								reciprocals)
{
	std::vector<char> identifier_seen(seq_db.size(), 0);
	unsigned cnt_line = 0;
	std::string line;
	std::vector<string> names;
	while (std::getline(f_seq_names, line))
	{
		++cnt_line;
		if (!line.length())
			continue;

		// if names are already in sequence file, must be indexed
		unsigned index;
		if (!seq_db.get_index_of_name(line, index))
		{
			std_print_error(CERR, "The sequence for \"" + line +
									"\" on line #" + ultoa(cnt_line) +
									" could not be found");
			continue;
		}

		// make sure each identifier only specified once
		if (identifier_seen[index])
		{
			std_print_error(CERR, "The identifier \"" + line +
									"\" has been specified already");
			continue;
		}
		identifier_seen[index] = 1;
		
		// all against this identifier
		if (reciprocals)
		{
			for (unsigned i = 0; i < seq_db.size(); ++i)
				if (i != index)
					sequence_pairs.push_back(make_pair(index, i));
		}		
		else
			for (unsigned i = 0; i < index; ++i)
				sequence_pairs.push_back(make_pair(index, i));
	}
	return true;
}



//________________________________________________________________________________________

//		parse_blast_pairs_from_file

//________________________________________________________________________________________
struct do_parse_pairs
{
	const t_collect_sequences&			seq_db;
	vector<pair<unsigned, unsigned> >& 	sequence_pairs;
	ostream&							CERR;
	string								comment_letters;
	do_parse_pairs(	const t_collect_sequences&			seq_db_,
					vector<pair<unsigned, unsigned> >& 	sequence_pairs_,
					ostream&							CERR_,
					const string&						comment_letters_):
			seq_db(seq_db_), sequence_pairs(sequence_pairs_), CERR(CERR_),
			comment_letters(comment_letters_){}
	void print_error(unsigned cnt_line, const string& prot_name) const
	{
		std_print_error(CERR, "The sequence for \"" + prot_name +
								"\" on line #" + ultoa(cnt_line) +
								" could not be found");
	}
	void operator()(istream& is, unsigned cnt_line)
	{
		// get data for one protein pair
		string prot_name1, prot_name2;

		// comment if false
		string ignore_comment_line;
		if (!seq_pairs::parse_ids(is, prot_name1, prot_name2, ignore_comment_line, comment_letters, seq_pairs::ePARSE_WHOLE_LINE))
			return;

		// if names are already in sequence file, must be indexed
		unsigned index1, index2;
		if (!seq_db.get_index_of_name(prot_name1,  index1))
		{
			print_error(cnt_line, prot_name1);
			return;
		}
		if (!seq_db.get_index_of_name(prot_name2,  index2))
		{
			print_error(cnt_line, prot_name1);
			return;
		}

		// found both names in sequence files
		sequence_pairs.push_back(make_pair(index1, index2));
	}
};





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



int main(int argc, char* argv[])
{
	debug_func_cerr( 8 );

	t_program_args args;

	try
	{

	// don't use C streams
	std::ios::sync_with_stdio(false);

	// get command line arguments
	if (!process_command_line_options(argc, argv, args))
		return 1;


	//
	// read sequences
	//
	t_collect_sequences seq_db(args.align_mode == eALIGN_PAIRWISE ? eSEQ_DUPLICATES : eSEQ_NODUPLICATES);
	{
		args.VERBOSE() << "\tReading sequences from STDIN...\n";
		t_progress_indicator dots(args.VERBOSE(), 10000);
		t_collect_sequences_report_errors_to_strm report_errors(args.CERR());
		seq_db.init_from_stream(std::cin, args.regex, dots, report_errors);
	}
	
	args.VERBOSE() << "\tDone...\n";
	args.VERBOSE() << "\t" << seq_db.size() << " sequences read.\n";


	vector<pair<unsigned, unsigned> > 	sequence_pairs;
	//
	// read sequences pairs
	//
	if (args.align_mode == eALIGN_SPECIFIED)
	{
		assert(args.file_pairs.length());
		args.VERBOSE() << "\tReading sequence pairs...\n";

		do_parse_pairs parse_pairs(seq_db,  sequence_pairs, args.CERR(), 
									args.comment_letters);
		seq_pairs::read_data(parse_pairs, args.file_pairs, args.VERBOSE());
	}

	//
	// all vs all
	//
	else if (args.align_mode == eALIGN_ALLVSALL)
	{
		args.VERBOSE() << "\tAll vs all comparisons...\n";
		unsigned sz = seq_db.size();
		if (args.reciprocals)
			for (unsigned i = 0; i < sz; ++i)
				for (unsigned j = 0; j < sz; ++j)
				{
					sequence_pairs.push_back(make_pair(i, j));
				}
		else if (args.self_self)
		{
			for (unsigned i = 0; i < sz; ++i)
				for (unsigned j = i; j < sz; ++j)
				{
					sequence_pairs.push_back(make_pair(i, j));
				}
		}
		else
			for (unsigned i = 0; i < sz - 1; ++i)
				for (unsigned j = i + 1; j < sz; ++j)
					sequence_pairs.push_back(make_pair(i, j));
	}
	else if (args.align_mode == eALIGN_IDENTIFIERS)
	{
		assert(args.file_identifiers.length());
		args.VERBOSE() << "\tReading identifiers to be compared against all other sequences...\n";
		ifstream f_seq_names(args.file_identifiers.c_str());
		if (!f_seq_names)
		{
			std_print_error(args.CERR(), "Could not open " + args.file_identifiers +
							" to retrieve identifiers");
			return 1;
		}
		read_seqid_identifiers(f_seq_names, seq_db, sequence_pairs, args.CERR(), args.reciprocals);
	}
	//
	// sequence pairs not specified: successive pairs
	//
	else
		//args.align_mode == eALIGN_PAIRWISE
	{
		args.VERBOSE() << "\tComparing successive pairs...\n";
		unsigned sz = seq_db.size();
		for (unsigned i = 0; i < sz - 1; i += 2)
			sequence_pairs.push_back(make_pair(i, i+1));
	}

	//
	//	blast options
	//
	t_portable_timer timer;
	t_blast_options blast_options(args.seg, args.cutoff, args.matrix, args.db_size);


	args.VERBOSE() << "\tAligning " << sequence_pairs.size() << " sequence pairs ...\n";


	if (args.results_header)
		cout	<< seq_pairs		::column_names() << "\t"
				<< t_align_results	::column_names() << "\n"
				<< seq_pairs		::column_name_underscores() << "\t"
				<< t_align_results	::column_name_underscores() << "\n";

	unsigned cnt_results = 0;
	{
		t_progress_indicator dots(args.VERBOSE(), std::max(20U, sequence_pairs.size() / 250),
																	sequence_pairs.size());
		dots.use_timer();

		//________________________________________________________________________________
		// 		align all genes in pairs
		//________________________________________________________________________________
		for (unsigned i = 0; i < sequence_pairs.size(); ++i)
		{
			unsigned seq1 = sequence_pairs[i].first;
			unsigned seq2 = sequence_pairs[i].second;
			if (align_sequences(dots, blast_options, args.lencutoff, seq_db, seq1, seq2, args))
				cnt_results++;
		}
	}

	if (cnt_results == sequence_pairs.size())
	{
		args.VERBOSE() << "\tAll " << cnt_results << " sequence pairs were aligned.\n\n";
	}
	else
	{
		std_print_warning(args.VERBOSE(), "Only " +
							ultoa(cnt_results) + " sequence pairs out of " +
							ultoa(sequence_pairs.size()) + " could be aligned.");
	}



	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}
	catch(...)
	{
		std_print_error(args.CERR(), "Software bug caused crash.\n");
		return 1;
	}

	return 0;
}

